#include <TROOT.h>

void Z4F_CHIRALSEA(int imom, int ioa2) 
{
        char Oa2print[8]; whichOa2(ioa2, Oa2print);
	char rgiprint[16];

        double ZAPREV=0.80;
        int Nmass;
	char dir[1024];

 	 /////////////////////////////
	// INPUT FILES///////////
        ////////////////////////////////


         // basic input file
        char input_dat[1024];
        sprintf(input_dat,"../INPUT/input_sea.dat");
//	printf("%s\n",input_dat);
        FILE *inputdat=fopen(input_dat,"r");
        fscanf(inputdat,"%i",&L);
        fscanf(inputdat,"%i",&T);
        fscanf(inputdat,"%i",&Nclust);
        fscanf(inputdat,"%s",dir);
	fscanf(inputdat,"%s",rgiprint);

        // mass input files

          char mPOL_dat[1024];
          sprintf(mPOL_dat,"../INPUT/mPOL_sea.dat");
//	  printf("%s\n",mPOL_dat);
          FILE *mPOLdat=fopen(mPOL_dat,"r");
          fscanf(mPOLdat,"%i",&Nmass);
          double *mPOL=(double*)malloc(sizeof(double)*(Nmass));
          int *e=(int*)malloc(sizeof(int)*(Nmass));
          for (int im=0; im<Nmass; im++)
		fscanf(mPOLdat,"%d",&e[im]);
          for (int im=0; im<Nmass; im++)
                fscanf(mPOLdat,"%lf",&mPOL[im]);


   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  Equivalent_Momenta(SinP2);

 //////////////////////////////////
 // IO files
  ///////////////////////////////

  char file[1024];
  char dir_input[1024];
  char dir_output[1024];
  sprintf(dir_output,"%s/a/QAV/",dir);

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

  //////////////////////////////////
  // Read Zs
   /////////////////////////////

     double ****ZpPC=malloc(sizeof(double***)*(5));
     for (int iop1=0; iop1<5; iop1++) ZpPC[iop1]=malloc(sizeof(double**)*(5));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC[iop1][iop2]=malloc(sizeof(double*)*(Nmass));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) for (int imass=0; imass<Nmass; imass++) ZpPC[iop1][iop2][imass]=malloc(sizeof(double)*(Nclust+1));

/*     double ****DeltapPC=malloc(sizeof(double***)*(5));
     for (int iop1=0; iop1<5; iop1++) DeltapPC[iop1]=malloc(sizeof(double**)*(5));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DeltapPC[iop1][iop2]=malloc(sizeof(double*)*(Nmass));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) for (int imass=0; imass<Nmass; imass++) DeltapPC[iop1][iop2][imass]=malloc(sizeof(double)*(Nclust+1));
*/

     double **Zq=malloc(sizeof(double*)*(Nmass));
     for (int imass=0; imass<Nmass; imass++) Zq[imass]=malloc(sizeof(double)*(Nclust+1));

     double *aux_clust=malloc(sizeof(double)*(Nclust+1));

     for (int imass=0; imass<Nmass; imass++)
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
     {
                sprintf(file,"ZpPC%d%d%s_%s",iop1+1,iop2+1,rgiprint,Oa2print);
		sprintf(dir_input,"%s/%da/QAV/",dir,e[imass]);
                read_input_onemomenta(dir_input,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) ZpPC[iop1][iop2][imass][iclust]=aux_clust[iclust];
     }
/*
     for (int imass=0; imass<Nmass; imass++)
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
     {
                sprintf(file,"DeltapPC%d%d_%s",iop1+1,iop2+1,Oa2print);
                sprintf(dir_input,"%s/%da/QAV/",dir,e[imass]);
                read_input_onemomenta(dir_input,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) DeltapPC[iop1][iop2][imass][iclust]=aux_clust[iclust];
     }
*/

    for (int imass=0; imass<Nmass; imass++)
    {
    sprintf(file,"Zq%s_%s",rgiprint,Oa2print);
    sprintf(dir_input,"%s/%da/QAV/",dir,e[imass]);
    read_input_onemomenta(dir_input,file,aux_clust,imom);
    for (int iclust=0; iclust<Nclust+1; iclust++) Zq[imass][iclust]=aux_clust[iclust];
    }

      ///////////////////////////////////////
      // CHIRAL FIT FOR EACH MOMENTUM 
     ///////////////////////////////////////

      double *ZAUX=malloc(sizeof(double)*(Nmass));
      double **ZAUX_CLUST=malloc(sizeof(double*)*(Nmass));
      for (int imass=0; imass<Nmass; imass++) ZAUX_CLUST[imass]=malloc(sizeof(double)*(Nclust+1));

      double *MAUX=malloc(sizeof(double)*(Nmass));

 
      double *ZqCHI=malloc(sizeof(double)*(Nclust+1));
      double *ZqLIN=malloc(sizeof(double)*(Nclust+1));
      double *ZqCHI2=malloc(sizeof(double)*(Nclust+1));

	

    // Chiral fit for Zq
    for (int iclust=0; iclust<Nclust+1; iclust++)
     {
     for (int imass=0; imass<Nmass; imass++) 
     {
                MAUX[imass]=mPOL[imass]*mPOL[imass];
                ZAUX[imass]=Zq[imass][iclust];
                for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=Zq[imass][iclust_loop];
     }
     lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,ZqCHI[iclust],ZqLIN[iclust],ZqCHI2[iclust]); 
     sprintf(file,"Zq%s_%s",rgiprint,Oa2print); write_output_sea_chiral_extrapolation(dir_output,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,ZqCHI,ZqLIN,ZqCHI2);

     }

      double ***ZCHI=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  ZCHI[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZCHI[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));      

      double ***ZLIN=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  ZLIN[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZLIN[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

      double ***ZCHI2=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  ZCHI2[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZCHI2[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

	

    // chiral fit of Z 


    for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
    for (int iclust=0; iclust<Nclust+1; iclust++)
     {
     for (int imass=0; imass<Nmass; imass++)
     {
                MAUX[imass]=mPOL[imass]*mPOL[imass];
                ZAUX[imass]=ZpPC[iop1][iop2][imass][iclust];
                for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=ZpPC[iop1][iop2][imass][iclust_loop];
     }
     lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,ZCHI[iop1][iop2][iclust],ZLIN[iop1][iop2][iclust],ZCHI2[iop1][iop2][iclust]);
     sprintf(file,"ZpPC%d%d%s_%s",iop1+1,iop2+1,rgiprint,Oa2print); write_output_sea_chiral_extrapolation(dir_output,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,ZCHI[iop1][iop2],ZLIN[iop1][iop2],ZCHI2[iop1][iop2]);
     }

/*
      double ***DeltapCHI=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  DeltapCHI[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DeltapCHI[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

      double ***DeltapLIN=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  DeltapLIN[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DeltapLIN[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

      double ***DeltapCHI2=malloc(sizeof(double**)*(5));
      for (int iop1=0; iop1<5; iop1++)  DeltapCHI2[iop1]=malloc(sizeof(double*)*(5));
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DeltapCHI2[iop1][iop2]=malloc(sizeof(double)*(Nclust+1));



    // chiral fit of Deltap 


    for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
    for (int iclust=0; iclust<Nclust+1; iclust++)
     {
     for (int imass=0; imass<Nmass; imass++)
     {
                MAUX[imass]=mPOL[imass]*mPOL[imass];
                ZAUX[imass]=DeltapPC[iop1][iop2][imass][iclust];
                for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=DeltapPC[iop1][iop2][imass][iclust_loop];
     }
     lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,DeltapCHI[iop1][iop2][iclust],DeltapLIN[iop1][iop2][iclust],DeltapCHI2[iop1][iop2][iclust]);
     sprintf(file,"DeltapPC%d%d%s_%s",iop1+1,iop2+1,rgiprint,Oa2print); write_output_sea_chiral_extrapolation(dir_output,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,DeltapCHI[iop1][iop2],DeltapLIN[iop1][iop2],DeltapCHI2[iop1][iop2]);

     }
*/
//********************************************************************************************************************///
///***************************************************************************************************************///////

    // write output
    printf("Writing output in %s/p%03d/(dat)(jk)\n",dir_output,imom);


   // Zq

   sprintf(file,"Zq%s_%s",rgiprint,Oa2print); 
   write_output_onemomenta_real(dir_output,file,ZqCHI,SinP2_eq,imom);
 
  //Z_ij
   for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
   {
       	sprintf(file,"ZpPC%d%d%s_%s",iop1+1,iop2+1,rgiprint,Oa2print);
	write_output_onemomenta_real(dir_output,file,ZCHI[iop1][iop2],SinP2_eq,imom);
   }
/* 
   for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
   {
        sprintf(file,"DeltapPC%d%d%s_%s",iop1+1,iop2+1,rgiprint,Oa2print);
        write_output_onemomenta_real(dir_output,file,DeltapCHI[iop1][iop2],SinP2_eq,imom);
   }
*/

    delete Zq; delete ZpPC; 
    delete aux_clust;
    delete ZAUX; delete ZAUX_CLUST; delete MAUX;
    delete ZqCHI; delete ZqLIN; delete ZqCHI2;
    delete ZCHI; delete ZLIN; delete ZCHI2; 
//    delete DeltapCHI; delete DeltapLIN; delete DeltapCHI2;


}
